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

Last change on this file since 6501 was 5165, checked in by ole, 16 years ago

Comments for bad_mesh

File size: 4.2 KB
Line 
1""" This Script shows how a bad mesh can be generated causing
2excessive values for stage and momentum.
3
4"""
5
6#------------------------------------------------------------------------------
7# Import necessary modules
8#------------------------------------------------------------------------------
9
10# Standard modules
11import os
12import time
13import sys
14
15# Related major packages
16from anuga.shallow_water import Domain
17from anuga.shallow_water import Dirichlet_boundary
18from anuga.pmesh.mesh_interface import create_mesh_from_regions
19from anuga.utilities.polygon import read_polygon, plot_polygons
20
21
22
23#------------------------------------------------------------------------------
24# Filenames
25#------------------------------------------------------------------------------
26
27outname ='bad_mesh'
28meshname ='bad_mesh.msh'
29
30
31#------------------------------------------------------------------------------
32# Create the triangular mesh based on overall clipping polygon with a tagged
33# boundary and interior regions defined in project.py along with
34# resolutions (maximal area of per triangle) for each polygon
35#------------------------------------------------------------------------------
36
37W=304180.0
38S=6185270.0
39E=307650.0
40N=6189040.0
41
42bounding_polygon = [[W, S], [E, S], [E, N], [W, N]]
43
44# Below is assigning the name for each polygon (roughness values added below)
45Refine_Read_Polygon1 = read_polygon('Catchment_Bdy.csv')
46Refine_Read_Polygon2 = read_polygon('OuterHarbour.csv')
47
48# Bad
49# interior_regions = [[Refine_Read_Polygon1, 10],[Refine_Read_Polygon2, 50]]
50
51# Bad (simpler) - Polygon 2 is warped
52interior_regions = [[Refine_Read_Polygon2, 50]]
53
54# OK
55#interior_regions = [[Refine_Read_Polygon1, 10]]
56
57# Output plot of polygons to file
58#plot_polygons([bounding_polygon,
59#               Refine_Read_Polygon1,
60#               Refine_Read_Polygon2],
61#               figname='milevsky')
62
63
64
65create_mesh_from_regions(bounding_polygon,
66                         boundary_tags={'south': [0],
67                                        'east': [1],
68                                        'north': [2],
69                                        'west': [3]},
70                         maximum_triangle_area=10000, #background mesh size
71                         interior_regions=interior_regions,
72                         filename=meshname,
73                         use_cache=False,
74                         verbose=True)
75
76
77#------------------------------------------------------------------------------
78# Setup computational domain
79#------------------------------------------------------------------------------
80
81domain = Domain(meshname, use_cache=False, verbose=True)
82domain.set_name(outname) 
83print domain.statistics()
84
85#------------------------------------------------------------------------------
86# Setup initial conditions
87#------------------------------------------------------------------------------
88# Setting up initial conditions for a Flood model includes:
89# Setting A Base topography over the model domain
90# Setting the variation of Roughness over the terrain, which may be
91# varied dependent on depth say
92# Setting initial water levels at BOundaries and in Lakes and Water Bodies.
93
94
95tide = 1.5 # 100yr tidal level as per WCC drainage design code
96base_friction = 0.015
97
98domain.set_quantity('friction', base_friction)
99domain.set_quantity('stage', tide)
100domain.set_quantity('elevation', 0)
101
102
103#------------------------------------------------------------------------------
104# Setup boundary conditions
105#------------------------------------------------------------------------------
106
107print 'Available boundary tags', domain.get_boundary_tags()
108
109Bd = Dirichlet_boundary([tide,0,0])
110
111# Boundary conditions
112domain.set_boundary({'west': Bd,
113                     'south': Bd,
114                     'north': Bd,
115                     'east': Bd})
116   
117
118#------------------------------------------------------------------------------
119# Evolve system through time
120#------------------------------------------------------------------------------
121
122import time
123t0 = time.time()
124
125for t in domain.evolve(yieldstep = 1, finaltime = 2): 
126    print domain.timestepping_statistics()
127
128print 'Finished'
129   
Note: See TracBrowser for help on using the repository browser.