source: anuga_work/debug/bad_mesh/bad_mesh.py @ 7244

""" This Script shows how a bad mesh can be generated causing
excessive values for stage and momentum.

"""

#------------------------------------------------------------------------------
# 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 Dirichlet_boundary
from anuga.pmesh.mesh_interface import create_mesh_from_regions
from anuga.utilities.polygon import read_polygon, plot_polygons



#------------------------------------------------------------------------------
# Filenames
#------------------------------------------------------------------------------

outname ='bad_mesh'
meshname ='bad_mesh.msh'


#------------------------------------------------------------------------------
# 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
#------------------------------------------------------------------------------

W=304180.0
S=6185270.0
E=307650.0
N=6189040.0

bounding_polygon = [[W, S], [E, S], [E, N], [W, N]]

# Below is assigning the name for each polygon (roughness values added below)
Refine_Read_Polygon1 = read_polygon('Catchment_Bdy.csv')
Refine_Read_Polygon2 = read_polygon('OuterHarbour.csv')

# Bad
# interior_regions = [[Refine_Read_Polygon1, 10],[Refine_Read_Polygon2, 50]]

# Bad (simpler) - Polygon 2 is warped
interior_regions = [[Refine_Read_Polygon2, 50]]

# OK 
#interior_regions = [[Refine_Read_Polygon1, 10]]

# Output plot of polygons to file
#plot_polygons([bounding_polygon,
#               Refine_Read_Polygon1,
#               Refine_Read_Polygon2],
#               figname='milevsky')



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


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

domain = Domain(meshname, use_cache=False, verbose=True)
domain.set_name(outname) 
print domain.statistics()

#------------------------------------------------------------------------------
# Setup initial conditions
#------------------------------------------------------------------------------
# Setting up initial conditions for a Flood model includes:
# Setting A Base topography over the model domain
# Setting the variation of Roughness over the terrain, which may be
# varied dependent on depth say
# Setting initial water levels at BOundaries and in Lakes and Water Bodies.


tide = 1.5 # 100yr tidal level as per WCC drainage design code
base_friction = 0.015

domain.set_quantity('friction', base_friction)
domain.set_quantity('stage', tide)
domain.set_quantity('elevation', 0)


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

print 'Available boundary tags', domain.get_boundary_tags()

Bd = Dirichlet_boundary([tide,0,0])

# Boundary conditions
domain.set_boundary({'west': Bd,
                     'south': Bd,
                     'north': Bd,
                     'east': Bd})
    

#------------------------------------------------------------------------------
# Evolve system through time
#------------------------------------------------------------------------------

import time
t0 = time.time()

for t in domain.evolve(yieldstep = 1, finaltime = 2): 
    print domain.timestepping_statistics()

print 'Finished'