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

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

example of adding segments and vertices to a mesh file, from a points file.

File size: 4.9 KB
Line 
1""" This Script shows how exta points and segments can be added to a msh file.
2
3# WARNING #
4* this doesn't have tests assosiated with it.  For example that the
5resolutions stay as described in create_mesh_from_regions is not tested
6* Also, this code causes the mesh to be generated twice, so it is not
7an efficient way of doing this.
8
9If this is used a lot you may want to wrap the functionality into a function.
10"""
11
12#------------------------------------------------------------------------------
13# Import necessary modules
14#------------------------------------------------------------------------------
15
16# Standard modules
17import os
18import time
19import sys
20
21# Related major packages
22from anuga.shallow_water import Domain
23from anuga.shallow_water import Dirichlet_boundary
24from anuga.pmesh.mesh_interface import create_mesh_from_regions
25from anuga.pmesh.mesh import importMeshFromFile
26from anuga.utilities.polygon import read_polygon, plot_polygons
27
28from anuga.geospatial_data.geospatial_data import Geospatial_data
29
30
31
32#------------------------------------------------------------------------------
33# Filenames
34#------------------------------------------------------------------------------
35
36outname ='mesh'
37meshname ='mesh.msh'
38
39# Create a points file
40ptfile = 'points_from_a_file.csv'
41file = open(ptfile,"w")
42file.write("x,y\n\
43304200.0, 6185300.0\n\
44304300.0, 6185400.0\n\
45304400.0, 6185500.0\n")
46file.close()
47
48#------------------------------------------------------------------------------
49# Create the triangular mesh based on overall clipping polygon with a tagged
50# boundary and interior regions defined in project.py along with
51# resolutions (maximal area of per triangle) for each polygon
52#------------------------------------------------------------------------------
53
54W=304180.0
55S=6185270.0
56E=307650.0
57N=6189040.0
58
59bounding_polygon = [[W, S], [E, S], [E, N], [W, N]]
60
61
62
63create_mesh_from_regions(bounding_polygon,
64                         boundary_tags={'south': [0],
65                                        'east': [1],
66                                        'north': [2],
67                                        'west': [3]},
68                         maximum_triangle_area=10000, #background mesh size
69                         filename=meshname,
70                         use_cache=False,
71                         verbose=True)
72
73m = importMeshFromFile(meshname)
74
75#m.add_vertices([(W +500,S +500),(W +1000,S +1000) ])
76
77# load points from a points file
78#m.add_vertices(Geospatial_data(ptfile))
79
80p1 = (W +500,N - 500)
81p2 = (W +1000,N - 1000)
82p3 = (W +1000,N - 1500) 
83#m.add_points_and_segments([p1,p2,p3], [[0,1],[1,2]])
84
85# add points and segments, where the segments are automatically added.
86# This will work on the next release of ANUGA.
87# It doesn't work on the current sourceforge release.
88p1 = (E -1000,N - 1000)
89p2 = (E -1000,N - 1500) 
90#m.add_points_and_segments([p1,p2])
91
92# load points from a points file
93m.add_points_and_segments(Geospatial_data(ptfile))
94
95m.generate_mesh()
96m.export_mesh_file('end.msh')
97
98import sys; sys.exit() 
99#------------------------------------------------------------------------------
100# Setup computational domain
101#------------------------------------------------------------------------------
102
103domain = Domain(meshname, use_cache=False, verbose=True)
104domain.set_name(outname) 
105print domain.statistics()
106
107#------------------------------------------------------------------------------
108# Setup initial conditions
109#------------------------------------------------------------------------------
110# Setting up initial conditions for a Flood model includes:
111# Setting A Base topography over the model domain
112# Setting the variation of Roughness over the terrain, which may be
113# varied dependent on depth say
114# Setting initial water levels at BOundaries and in Lakes and Water Bodies.
115
116
117tide = 1.5 # 100yr tidal level as per WCC drainage design code
118base_friction = 0.015
119
120domain.set_quantity('friction', base_friction)
121domain.set_quantity('stage', tide)
122domain.set_quantity('elevation', 0)
123
124
125#------------------------------------------------------------------------------
126# Setup boundary conditions
127#------------------------------------------------------------------------------
128
129print 'Available boundary tags', domain.get_boundary_tags()
130
131Bd = Dirichlet_boundary([tide,0,0])
132
133# Boundary conditions
134domain.set_boundary({'west': Bd,
135                     'south': Bd,
136                     'north': Bd,
137                     'east': Bd})
138   
139
140#------------------------------------------------------------------------------
141# Evolve system through time
142#------------------------------------------------------------------------------
143
144import time
145t0 = time.time()
146
147for t in domain.evolve(yieldstep = 1, finaltime = 2): 
148    print domain.timestepping_statistics()
149
150print 'Finished'
151   
Note: See TracBrowser for help on using the repository browser.