Changeset 5569 for anuga_work/production/perth/run_perth.py

Timestamp:

Jul 25, 2008, 8:53:51 AM (15 years ago)

Author:

kristy

Message:

updated scripts for polyline boundary

File:

• anuga_work/production/perth/run_perth.py (modified) (6 diffs)

anuga_work/production/perth/run_perth.py

@@ -30 +30 @@
 from anuga.shallow_water import Field_boundary
 from Numeric import allclose
-from anuga.shallow_water.data_manager import export_grid
+from anuga.shallow_water.data_manager import export_grid, create_sts_boundary
 
 from anuga.pmesh.mesh_interface import create_mesh_from_regions
@@ -39 +39 @@
 from anuga.damage_modelling.inundation_damage import add_depth_and_momentum2csv, inundation_damage
 from anuga.fit_interpolate.benchmark_least_squares import mem_usage 
+from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
 
 # Application specific imports
@@ -70 +71 @@
     print "Processor Name:",get_processor_name()
 
-    # filenames
-#    meshes_dir_name = project.meshes_dir_name+'.msh'
-
-    # creates copy of code in output dir
-    print 'min triangles', project.trigs_min,
-    print 'Note: This is generally about 20% less than the final amount'
-
+    
+    #-----------------------------------------------------------------------
+    # Domain definitions
+    #-----------------------------------------------------------------------
+
+    # Read in boundary from ordered sts file
+    urs_bounding_polygon=create_sts_boundary(project.scenario_name)
+
+    # Reading the landward defined points, this incorporates the original clipping
+    # polygon minus the 100m contour
+    landward_bounding_polygon = read_polygon(project.polygons_dir+'landward_bounding_polygon.txt')
+
+    # Combine sts polyline with landward points
+    bounding_polygon = urs_bounding_polygon + landward_bounding_polygon
+    
+    # counting segments
+    N = len(urs_bounding_polygon)-1
+    boundary_tags={'back': [N+1,N+2,N+3], 'side': [N,N+4],'ocean': range(N)}
+
+    
     #--------------------------------------------------------------------------
     # Create the triangular mesh based on overall clipping polygon with a
@@ -90 +104 @@
         print 'start create mesh from regions'
 
-        create_mesh_from_regions(project.poly_all,
-                             boundary_tags=project.boundary_tags,
+        create_mesh_from_regions(bounding_polygon,
+                             boundary_tags=boundary_tags,
                              maximum_triangle_area=project.res_poly_all,
                              interior_regions=project.interior_regions,
@@ -168 +182 @@
     print 'domain id', id(domain)
    
-
+    boundary_urs_out=project.boundaries_dir_name
+    
     print 'Available boundary tags', domain.get_boundary_tags()
     Bf = File_boundary(boundary_urs_out+'.sts', 
@@ -174 +189 @@
                    use_cache=True,
                    verbose = True,
-                   boundary_polygon=project.bounding_polygon)
+                   boundary_polygon=bounding_polygon)
 
     Br = Reflective_boundary(domain)