source: anuga_work/production/sydney_2006/eq_test.py @ 4344

from anuga.shallow_water import Domain, Reflective_boundary, Dirichlet_boundary 
from eqf_v2 import earthquake_tsunami

res = 10000000.
intres = 1000000.
boundingpoly=[[-100000,-80000],[-100000,80000],[100000,80000],[100000,-80000]]
poly = [[-40000,-20000],[40000,-20000],[40000,30000],[-40000,30000]]
interior_regions = [[poly,intres]]

from anuga.pmesh.mesh_interface import create_mesh_from_regions
from caching import cache
meshname = 'another_test'+'.msh'

_ = cache(create_mesh_from_regions,
          boundingpoly,
           {'boundary_tags': {'e0': [0], 'e1': [1], 'e2': [2], 'e3': [3]},
           'maximum_triangle_area': res,
           'filename': meshname,
           'interior_regions': interior_regions},
           verbose = True)

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

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

tsunami_source = earthquake_tsunami(length=25.0,
                                    width=5.,
                                    strike=0.,
                                    depth=3500.0,
                                    dip=15.0,
                                    slip=10.,
                                    x0=0.0,
                                    y0=0.0,
                                    domain=domain,
                                    verbose=True)

stage0 = 1.0
#domain.set_quantity('stage', stage0)
#z1 = domain.get_quantity('stage')
#int1 = z1.get_integral()/400000.0

domain.set_quantity('stage', tsunami_source)
z2 = domain.get_quantity('stage')
int2 = z2.get_integral()/400000.0
#print 'integral before and after', int1, int2

print len(z2)
print dir(z2)
#X,Y,A,V = z2.get_vertex_values()
#print 'x', V

#print 'values', z.get_values()
#print 'vertex', z.get_vertex_values()
#import sys; sys.exit()

domain.set_quantity('elevation', -10.0, alpha = 0.1)

Br = Reflective_boundary(domain)
Bd = Dirichlet_boundary([0,0,0])
domain.set_boundary( {'e0': Br,  'e1': Br, 'e2': Br, 'e3': Br} )

import time

t0 = time.time()

for t in domain.evolve(yieldstep = 1, finaltime = 1): 
    domain.write_time()
    domain.write_boundary_statistics(tags = 'e2')