Changeset 3799

Timestamp:

Oct 16, 2006, 5:47:09 PM (18 years ago)

Author:

ole

Message:

Work on okushiri mesh generation to use new mesh_interface

Location:

anuga_validation/automated_validation_tests/okushiri_tank_validation

Files:

• create_okushiri.py (modified) (7 diffs)
• validate_okushiri.py (modified) (4 diffs)

anuga_validation/automated_validation_tests/okushiri_tank_validation/create_okushiri.py

@@ -2 +2 @@
 """
 
-#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-# Assume that the root AnuGA dir is included in your pythonpath
-#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 from Numeric import array, zeros, Float, allclose
 
 from anuga.pmesh.mesh import *
+from anuga.pmesh.mesh_interface import create_mesh_from_regions
 from anuga.coordinate_transforms.geo_reference import Geo_reference
 
@@ -73 +71 @@
 if __name__ == "__main__":
 
-    #Preparing points
-    #(Only when using original Benchmark_2_Bathymetry.txt file)
-    #from anuga.abstract_2d_finite_volumes.data_manager import xya2pts
-    #xya2pts(project.bathymetry_filename, verbose = True,
-    #        z_func = lambda z: -z)
-
-
-
-
-
     # Prepare time boundary
     prepare_timeboundary(project.boundary_filename)
 
 
-    # Create Mesh
-
-    #Basic geometry
-    
+    #--------------------------------------------------------------------------
+    # Create the triangular mesh based on overall clipping polygon with a
+    # tagged
+    # boundary and interior regions defined in project.py along with
+    # resolutions (maximal area of per triangle) for each polygon
+    #--------------------------------------------------------------------------
+
+    base_resolution = 1
+
+    # Basic geometry
     xleft   = 0
     xright  = 5.448
@@ -96 +90 @@
     ytop    = 3.402
 
-    #Outline 
+    # Outline 
     point_sw = [xleft, ybottom]
     point_se = [xright, ybottom]
@@ -102 +96 @@
     point_ne = [xright, ytop]
 
-    #Midway points (left)
+    # Midway points (left)
     point_mtop = [xleft + (xright-xleft)/3+1, ytop]
     point_mbottom = [xleft + (xright-xleft)/3+1, ybottom]
 
-
-    geo = Geo_reference(xllcorner = xleft,
-                        yllcorner = ybottom)
-    
-                        
-    print "***********************"
-    print "geo ref", geo
-    print "***********************"
-    
-    m = Mesh(geo_reference=geo)
-
-    #Boundary
-    dict = {}
-    dict['points'] = [point_se,   #se
-                      point_ne,
-                      point_mtop,
-                      point_nw,
-                      point_sw,
-                      point_mbottom]
-
-    
-    dict['segments'] = [[0,1], [1,2], [2,3], [3,4],
-                        [4,5], [5,0],  #The outer border
-                        [2,5]]         #Separator
-    
-    dict['segment_tags'] = ['wall',
-                            'wall',
-                            'wall',
-                            'wave',
-                            'wall',
-                            'wall',
-                            '']           #Interior
-
-        
-    m.addVertsSegs(dict)
-
-
+    #Localised refined area for gulleys
+    xl = 4.8
+    xr = 5.3
+    yb = 1.6
+    yt = 2.3
+    p0 = [xl, yb]
+    p1 = [xr, yb]
+    p2 = [xr, yt]
+    p3 = [xl, yt]
+    gulleys = [p0, p1, p2, p3]
+
+    #Island area
+    island_0 = [xleft + 2*(xright-xleft)/3+1.2, ytop-0.5]
+    island_1 = [xleft + 2*(xright-xleft)/3+0.5, ybottom + 2*(ytop-ybottom)/3]
+    island_2 = [xleft + (xright-xleft)/2+0.3, ybottom + 2*(ytop-ybottom)/3-0.3]
+    island_3 = [xleft + (xright-xleft)/2+0.3, ybottom + (ytop-ybottom)/3+0.3]
+    island_4 = [xleft + 2*(xright-xleft)/3+0.4, ybottom + (ytop-ybottom)/3-0.3]
+    island_5 = [xleft + 2*(xright-xleft)/3+1.2, ybottom+0.2]
+    island_6 = [xl-.01, yb]  #OK
+    island_7 = [xl-.01, yt]  #OK      
+
+    island = [island_0, island_1, island_2,
+              island_3, island_4, island_5,
+              island_6, island_7]
+
+
+    
+
+    bounding_polygon = [point_se,
+                        point_ne,
+                        point_mtop,
+                        point_nw,
+                        point_sw,
+                        point_mbottom]
+    
+
+    interior_regions = [[gulleys, 0.00002*base_resolution],
+                        [island, 0.0003*base_resolution]]
+    
+
+
+    
+
+    print 'start create mesh from regions'
+
+    meshname = project.mesh_filename + '.msh'
+    m = create_mesh_from_regions(bounding_polygon,
+                                 boundary_tags={'wall': [0, 1, 2, 4, 5],
+                                                'wave': [3]},
+                                 maximum_triangle_area=0.1,
+                                 interior_regions=interior_regions,
+                                 use_cache=True,
+                                 verbose=True)
+
+    m.generateMesh('pzq28.0za1000000a') #???
+    
+    import project
+    m.export_mesh_file(project.mesh_filename)
+
+    import sys; sys.exit() 
 
 
@@ -175 +192 @@
     dict['points'] = [island_0, island_1, island_2,
                       island_3, island_4, island_5,
-                      #p0, p3]                      
                       island_6, island_7]
 
@@ -190 +206 @@
 #
     
-    base_resolution = 1
+    base_resolution = 10
 
     ocean = m.addRegionEN(xleft+.1, ybottom+.1)
@@ -207 +223 @@
 
 
-    # From r 1709 11 August 2005
-    #base_resolution = 100
-    #
-    #ocean = m.addRegionEN(xleft+.1, ybottom+.1)
-    #ocean.setMaxArea(0.01*base_resolution)
-    #
-    #mid = m.addRegionEN(point_mbottom[0]+.1, ybottom+.1)
-    #mid.setMaxArea(0.001*base_resolution)
-
-    #inner = m.addRegionEN(island_3[0]+.1, island_3[1]+.1)
-    #inner.setMaxArea(0.0001*base_resolution)
-
-
-    #gulleys = m.addRegionEN(p0[0]+0.1, p0[1]+0.1)
-    #gulleys.setMaxArea(0.00001*base_resolution)        
-    
-    
     m.generateMesh('pzq28.0za1000000a')
 

anuga_validation/automated_validation_tests/okushiri_tank_validation/validate_okushiri.py

@@ -37 +37 @@
 #-------------------------
 #N = 50
-N = 50
-points, vertices, boundary = rectangular_cross(N, N/5*3,
-                                               len1=5.448, len2=3.402)
-domain = Domain(points, vertices, boundary)
-
-domain.set_name('okushiri_automated_validation')
+#N = 50
+#points, vertices, boundary = rectangular_cross(N, N/5*3,
+#                                               len1=5.448, len2=3.402)
+#domain = Domain(points, vertices, boundary)
+
+
+print 'Creating domain from', project.mesh_filename
+
+domain = Domain(project.mesh_filename, use_cache=True, verbose=True)
+
+
+domain.set_name('okushiri_automated_validation_varmesh')
 domain.set_default_order(2)
 domain.set_minimum_storable_height(0.001)
+domain.set_maximum_allowed_speed(0) # The default in August 2005
+
 domain.check_integrity()
 
@@ -67 +75 @@
                          domain,
                          verbose=True)
+
 Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
-domain.set_boundary({'left': Bts, 'right': Br, 'bottom': Br, 'top': Br})
+
+domain.set_boundary({'wave': Bts, 'wall': Br})
+
+
 
 
@@ -127 +139 @@
 t0 = time.time()
 
+w_max = 0
 for j, t in enumerate(domain.evolve(yieldstep = timestep,
                                     finaltime = finaltime)):
@@ -139 +152 @@
         #print 'difference', x, y,\
         #      predicted_gauge_values[i][j] - reference_gauge_values[i][j]
+
+
+    w = domain.get_maximum_inundation_elevation()
+    x, y = domain.get_maximum_inundation_location()
+    print '  Coastline elevation = %.2f at (x,y)=(%.2f, %.2f)' %(w, x, y)
+    print
+    
+    if w > w_max:
+        w_max = w
+        x_max = x
+        y_max = y
+
+
+print '**********************************************'
+print 'Max coastline elevation = %.2f at (%.2f, %.2f)' %(w_max, x_max, y_max)
+print '**********************************************'
+