source: anuga_work/development/cocos/find_alpha_cocos.py @ 4991

Last change on this file since 4991 was 4991, checked in by nick, 16 years ago

add find_alpha_cocos.py to help fix TRAC ticket 230

File size: 3.7 KB
Line 
1
2
3from anuga.geospatial_data.geospatial_data import *
4from os import sep, environ, getenv, getcwd,umask
5from anuga.pmesh.mesh_interface import create_mesh_from_regions
6from anuga.shallow_water import Domain
7from anuga.shallow_water import Reflective_boundary
8from anuga.shallow_water.data_manager import export_grid
9from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files
10from math import sqrt
11from time import localtime, strftime, gmtime
12from anuga.utilities.polygon import number_mesh_triangles
13
14
15umask(002)
16time = strftime('%Y%m%d_%H%M%S',gmtime())
17
18home = getenv('INUNDATIONHOME') + sep +'data'+sep
19#output_dir = home + 'anuga_validation'+sep+'cocos'+sep+'anuga'+sep+'outputs'+sep+time+sep
20output_dir = '/d/ehn/4/ehn/tsunami/data/bathymetry/Cocos/bathymetry'+sep+time+sep
21
22#copy_code_files(output_dir,__file__)
23
24#start_screen_catcher(output_dir,verbose=True)
25print'testing'
26
27data_dir='/d/ehn/4/ehn/tsunami/data/bathymetry/Cocos/bathymetry/'
28data_dir_name = data_dir + 'new_edited_cocos_final'
29
30data_dir_name_small=data_dir_name+'_small'
31data_dir_name_smaller=data_dir_name+'_smaller'
32data_dir_name_small_clip=data_dir_name+'_small_clip'
33data_dir_name_clip=data_dir_name+'_clip'
34#print'create Geospatial data2 objects from coast'
35#G = Geospatial_data(file_name = data_dir_name + '.txt')
36#print'start split'
37#G_small, G_other = G.split(factor=0.1, verbose=True)
38#print'start export'
39#G_small.export_points_file(data_dir_name_smaller + '.txt')
40
41mesh_file=data_dir+'cocos_alpha.msh'
42sww_file='test_final'
43
44poly1 = [[439850.00,8565314.80],[439850.00, 8783110.00],
45        [156494.00,8783110.00],[156494.00,8565314.80]]
46#        68450m2
47grid1_res=370.4
48
49
50poly2 = [[294652.80,8625096.92],[294652.80,8713178.04],
51         [224424.96,8713178.04],[224424.96,8625096.92]]
52#        2738m2
53grid2_res=74.08
54
55poly3 = [[274295.91,8657346.91],[274295.91,8665112.56],
56         [266767.44,8665112.56],[266767.44,8657346.91]]
57#        110m2
58grid3_res=14.82
59       
60poly4 = [[270900.24,8659120.05],[270900.24,8660170.84],
61         [269701.45,8660170.84],[269701.45,8659120.05]]
62#        4.3m2
63grid4_res=2.96
64
65res = 100
66#G = Geospatial_data(file_name = data_dir_name + '.txt')
67#G1=G.clip(poly1)
68#G1.export_points_file(data_dir_name_clip + '.txt')
69
70
71interior_polys = [[poly2,(grid2_res**2/2)*res],[poly3,(grid3_res**2/2)*res],[poly4,(grid4_res**2/2)*res]]
72
73
74trigs_min = number_mesh_triangles(interior_polys, poly1, (grid1_res**2/2)*res)
75
76create_mesh_from_regions(poly1,
77                         boundary_tags={'side': [0,1,2,3]},
78                         maximum_triangle_area=(grid1_res**2/2)*res,
79                         interior_regions=interior_polys,
80                         filename=mesh_file,
81                         use_cache=True,
82                         verbose=True)
83                         
84                         
85#value, alpha = find_optimal_smoothing_parameter(data_file=data_dir_name_clip + '.txt',
86value, alpha = find_optimal_smoothing_parameter(data_file=data_dir_name_small + '.txt', 
87#                                             alpha_list=[0.0001,0.001,0.01,0.1,1.0],
88                                             alpha_list=[0.01,0.05,0.1,0.5,1.0],
89 #                                            alpha_list=[0.1],
90                                             boundary_poly=poly1,
91                                             mesh_file=mesh_file,
92                                             plot_name=None,
93                                             split_factor=0.001,
94                                             cache=True,
95                                             verbose=True)
96
97                                             
Note: See TracBrowser for help on using the repository browser.