# Changeset 659

Ignore:
Timestamp:
Dec 2, 2004, 2:11:36 PM (19 years ago)
Message:

Added functionality for settimng quantities according to polygons

Location:
inundation/ga/storm_surge
Files:
7 edited

Unmodified
Removed
• ## inundation/ga/storm_surge/examples/beach.py

 r653 from pmesh2domain import pmesh_to_domain_instance from config import data_dir from util import read_polygon, Polygon_function ###################### def inside(p,polygon): return True def friction(x, y): """Read friction from polygon file and set accordingly """ from Numeric import array x = array(x).astype('f') y = array(y).astype('f') p0 = [[20,20], [40,20], [40,40], [20,40]] #Get polygon fid = open('beach.poly') lines = fid.readlines() fid.close() polygon = [] for line in lines: fields = line.split(',') polygon.append( [float(fields[0]), float(fields[1])] ) print polygon z = 0*x for i in range(len(x)): if inside( (x[i], y[i]), polygon): z[i] = 1 else: z[i] = 0 friction([0],[1]) domain.set_quantity('elevation', x_slope) domain.set_quantity('level', x_slope) domain.set_quantity('elevation', Polygon_function( [(p0, 50)]), 'centroids' ) domain.set_quantity('level', Polygon_function( [(p0, 50)] ), 'centroids' ) domain.set_quantity('friction', 0.07) print domain.get_extent() import sys; sys.exit() ######################
• ## inundation/ga/storm_surge/pyvolution/domain.py

 r648 self.set_quantity(key, quantity_dict[key], location='vertices') def set_quantity(self, name, X, location='vertices', indexes = None): """Set values for named quantity
• ## inundation/ga/storm_surge/pyvolution/flatbed.py

 r505 from mesh_factory import rectangular from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\ Constant_height, wind_Stress Constant_height from Numeric import array from util import Polygon_function, read_polygon #Create basic mesh N = 150 points, vertices, boundary = rectangular(N, N, 1000, 1000) N = 20 points, vertices, boundary = rectangular(N, N, 100, 100) #Create shallow water domain domain = Domain(points, vertices, boundary) domain.smooth = True domain.visualise = False domain.store = True domain.default_order=1 domain.set_name('polygons') domain.default_order=2 #Set driving forces inflow_level = 10.0 domain.set_quantity('friction', manning) #Define polynomials p0 = [[20,27], [30,25], [40,40], [20,40]] p1 = [[80,19], [90,20], [85,50], [80,55], [75,58], [70,60], [60,24]] p2 = read_polygon('testpoly.txt') #Set elevation domain.set_quantity('elevation', Polygon_function([(p0, 23), (p1,10), (p2,15)])) ###################### #Evolution for t in domain.evolve(yieldstep = 10, finaltime = 500): for t in domain.evolve(yieldstep = 1, finaltime = 1): domain.write_time()
• ## inundation/ga/storm_surge/pyvolution/quantity.py

 r593 #Intialise centroid and edge_values self.interpolate() if location == 'centroids': #Extrapolate 1st order - to capture notion of area being specified self.extrapolate_first_order() def get_values(self, location='vertices', indexes = None): def extrapolate_first_order(self): """Extrapolate conserved quantities from centroid to vertices for each volume using first order scheme. """ qc = self.centroid_values qv = self.vertex_values for i in range(3): qv[:,i] = qc boundary values, storage and method for updating, and methods for extrapolation from centropid to vertices inluding methods for (second order) extrapolation from centroid to vertices inluding gradients and limiters """ def extrapolate_first_order(self): """Extrapolate conserved quantities from centroid to vertices for each volume using first order scheme. """ qc = self.centroid_values qv = self.vertex_values for i in range(3): qv[:,i] = qc def extrapolate_second_order(self): #Call correct module function
• ## inundation/ga/storm_surge/pyvolution/shallow_water.py

 r623 bed = self.quantities['elevation'] msg = 'All water levels must be greater than the bed elevation' assert alltrue( greater_equal( level.vertex_values, bed.vertex_values )), msg assert alltrue( greater_equal( level.edge_values, bed.edge_values )), msg assert alltrue( greater_equal( level.centroid_values, bed.centroid_values )), msg #msg = 'All water levels must be greater than the bed elevation' #assert alltrue( greater_equal( #    level.vertex_values, bed.vertex_values )), msg # #assert alltrue( greater_equal( #    level.edge_values, bed.edge_values )), msg # #assert alltrue( greater_equal( #    level.centroid_values, bed.centroid_values )), msg
• ## inundation/ga/storm_surge/pyvolution/test_util.py

 r656 def test_point_on_line(self): #Endpoints first assert point_on_line( (0, 0), [(0,0), (1,0)] ) assert not point_on_line( (0, 0.5), [(0,0), (1,1)] ) #From real example that failed assert not point_on_line( (40,50), [(40,20), (40,40)] ) #From real example that failed assert not point_on_line( (40,19), [(40,20), (40,40)] ) def test_inside_polygon(self): #Simplest case: Polygon is the unit square polygon = [[0,0], [1,0], [1,1], [0,1]] assert not inside_polygon( (0.5, 0.), polygon, closed=False) assert not inside_polygon( (1., 0.5), polygon, closed=False) #From real example (that failed) polygon = [[20,20], [40,20], [40,40], [20,40]] points = [ [40, 50] ] res = inside_polygon(points, polygon) assert len(res) == 0 polygon = [[20,20], [40,20], [40,40], [20,40]] points = [ [25, 25], [30, 20], [40, 50], [90, 20], [40, 90] ] res = inside_polygon(points, polygon) assert len(res) == 2 assert allclose(res, [0,1]) #More convoluted and non convex polygon #assert not inside_polygon( (0.5, 0.5), polygon ) points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
• ## inundation/ga/storm_surge/pyvolution/util.py

 r656 y1 = line[1][1] a = array([x - x0, y - y0]) a = array([x - x0, y - y0]) a_normal = array([a[1], -a[0]]) b = array([x1 - x0, y1 - y0]) len_a = sqrt(sum(a**2)) len_b = sqrt(sum(b**2)) if len_a == 0 or len_a == len_b: return True if dot(a_normal, b) == 0: #Point is somewhere on the infinite extension of the line len_a = sqrt(sum(a**2)) len_b = sqrt(sum(b**2)) if dot(a, b) >= 0 and len_a <= len_b: return True else: return False else: x = dot(a, b)/len_b return allclose(x, len_a) return False def inside_polygon(point, polygon, closed = True): """Determine whether points are inside or outside a polygon py = polygon[:,1] #Begin algorithm indices = [] x = points[k, 0] y = points[k, 1] inside = False for polygon, value in self.regions: indices = inside_polygon(points, polygon) for i in indices: z[i] = value #from sys import exit; exit() return z def read_polygon(filename): """Read points assumed to form a polygon There must be exactly two numbers in each line """ #Get polygon fid = open(filename) lines = fid.readlines() fid.close() polygon = [] for line in lines: fields = line.split(',') polygon.append( [float(fields[0]), float(fields[1])] ) return polygon class File_function:
Note: See TracChangeset for help on using the changeset viewer.