source: anuga_work/debug/will_power_example/make_mesh.py @ 5764

Last change on this file since 5764 was 5764, checked in by duncan, 16 years ago

An example of adding points and segments to a mesh.

File size: 4.5 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.pmesh.mesh import importMeshFromFile
20from anuga.utilities.polygon import read_polygon, plot_polygons
21
22from anuga.geospatial_data.geospatial_data import Geospatial_data
23
24
25
26#------------------------------------------------------------------------------
27# Filenames
28#------------------------------------------------------------------------------
29
30outname ='mesh'
31meshname ='mesh.msh'
32
33# Create a points file
34ptfile = 'points_from_a_file.csv'
35file = open(ptfile,"w")
36file.write("x,y\n\
37304200.0, 6185300.0\n\
38304250.0, 6185350.0\n\
39304300.0, 6185400.0\n")
40file.close()
41
42#------------------------------------------------------------------------------
43# Create the triangular mesh based on overall clipping polygon with a tagged
44# boundary and interior regions defined in project.py along with
45# resolutions (maximal area of per triangle) for each polygon
46#------------------------------------------------------------------------------
47
48W=304180.0
49S=6185270.0
50E=307650.0
51N=6189040.0
52
53bounding_polygon = [[W, S], [E, S], [E, N], [W, N]]
54
55
56
57create_mesh_from_regions(bounding_polygon,
58                         boundary_tags={'south': [0],
59                                        'east': [1],
60                                        'north': [2],
61                                        'west': [3]},
62                         maximum_triangle_area=10000, #background mesh size
63                         filename=meshname,
64                         use_cache=False,
65                         verbose=True)
66
67m = importMeshFromFile(meshname)
68
69m.add_vertices([(W +500,S +500),(W +1000,S +1000) ])
70
71# load points from a points file
72m.add_vertices(Geospatial_data(ptfile))
73
74p1 = (W +500,N - 500)
75p2 = (W +1000,N - 1000)
76p3 = (W +1000,N - 1500) 
77m.add_points_and_segments([p1,p2,p3], [[0,1],[1,2]])
78
79# add points and segments, where the segments are automatically added.
80# This will work on the next release of ANUGA.
81# It doesn't work on the current sourceforge release.
82p1 = (E -1000,N - 1000)
83p2 = (E -1000,N - 1500) 
84#m.add_points_and_segments([p1,p2])
85
86m.generate_mesh()
87m.export_mesh_file('end.msh')
88
89import sys; sys.exit() 
90#------------------------------------------------------------------------------
91# Setup computational domain
92#------------------------------------------------------------------------------
93
94domain = Domain(meshname, use_cache=False, verbose=True)
95domain.set_name(outname) 
96print domain.statistics()
97
98#------------------------------------------------------------------------------
99# Setup initial conditions
100#------------------------------------------------------------------------------
101# Setting up initial conditions for a Flood model includes:
102# Setting A Base topography over the model domain
103# Setting the variation of Roughness over the terrain, which may be
104# varied dependent on depth say
105# Setting initial water levels at BOundaries and in Lakes and Water Bodies.
106
107
108tide = 1.5 # 100yr tidal level as per WCC drainage design code
109base_friction = 0.015
110
111domain.set_quantity('friction', base_friction)
112domain.set_quantity('stage', tide)
113domain.set_quantity('elevation', 0)
114
115
116#------------------------------------------------------------------------------
117# Setup boundary conditions
118#------------------------------------------------------------------------------
119
120print 'Available boundary tags', domain.get_boundary_tags()
121
122Bd = Dirichlet_boundary([tide,0,0])
123
124# Boundary conditions
125domain.set_boundary({'west': Bd,
126                     'south': Bd,
127                     'north': Bd,
128                     'east': Bd})
129   
130
131#------------------------------------------------------------------------------
132# Evolve system through time
133#------------------------------------------------------------------------------
134
135import time
136t0 = time.time()
137
138for t in domain.evolve(yieldstep = 1, finaltime = 2): 
139    print domain.timestepping_statistics()
140
141print 'Finished'
142   
Note: See TracBrowser for help on using the repository browser.